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ABSTRACT 

A nonlinear elastic force-displacement relationship is used to calculate the 
transient impact force and local deformation at the point of contact between impactor and 
target. The nonlinear analysis and transfer function capabilities of NASTRAN are used to 
define a finite element model that behaves globally linearly elastic, and locally nonlinear 
elastic to model the local contact behavior. Results are presented for two different 
structures: a uniform cylindrical rod impacted longitudinally; and an orthotropic plate 
impacted transversely. Calculated impact force and transient structural response of the 
targets are shown to compare well with results measured in experimental tests. 


INTRODUCTION 

Aerospace structures are subjected to impact loading from a variety of sources, 
including dropped tools, runway debris, and munitions. In advanced composite materials, 
impact loading can cause significant internal structural damage, and a resulting loss of 
stiffness and strength. Therefore, the development of an accurate means of calculating 
structural response due to impact loading is of critical importance in the analysis and 
design of advanced aerospace structures. In this paper, a computational technique is 
developed to predict the dynamic response of a structure to a low velocity elastic impact. 

Current work in this area has led to a variety of methods for calculating the 
response of a structure subject to impact loading. Lai [1] used 3 equivalent springs in series 
to model impact of a composite plate. The springs represented the indentation, and 
flexural and shear stiffnesses. The calculated displacement results compared well with his 
test measurements for low velocity impact. Schonberg, et al. [2] developed a dosed form 
solution for the transverse impact of beams and plates by superimposing a static layer 
solution with elementary plate theory. His results showed that the impact force, maximum 
displacement and duration of contact were directly proportional to the target’s mass, 
although the analysis could not be readily applied to inhomogeneous materials. Graves [3] 
wrote a FORTRAN computer program to model the impact of composite plates. He 
predicted damage contours in the target using a maximum shear stress criterion. Lee, et al. 
[4] developed a specialty finite element code and used a triangular pulse load to simulate 
the impact force applied to a composite laminate. They showed that the higher modes of 
the target’s frequency response are more pronounced near the point of contact with the 
impactor. Sun [5] proposed a modified nonlinear elastic contact law for orthotropic 
materials that accounted for permanent indentation during unloading. The new contact 
law was used in a finite element program to calculate the dynamic response of composite 
beams under transverse impact loading. Tan [6] extended this approach to laminated 
plates. He used a nine— node isoparametric plate element to model the target, and his 
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calculated strain and impact force histories showed good agreement with test data. Sun [7] 
then included the effects of in-plane prestress on the plate. The calculated transient strain 
response in the target laminate showed good agreement with their experimental 
measurements. Aggour [8] extended this work to include the effects of transverse shear and 
rotary inertia. Strain response from their analysis compared well with experimental 
measurements obtained by Takeda [9-10]. Takamatsu, et al [11] authored a three 
dimensional finite element code to predict impact damage in composite laminates. A 
maximum stress failure criterion, as developed by Tsai [12], was used as a basis for 
progressively altering the stiffness of the composite target. Humphreys [13] developed a 
finite element code that also included the effects of material damping, shear deformation, 
and delamination. He used a stress-based failure criterion. Wu and Springer [14] 
developed a 3 dimensional transient finite element code that predicts delamination damage 
based on a novel application of dimensional analysis. Calculated delamination contours 
were shown to compare favorably with ultrasonic c— scans of impacted test specimens. 

The most common approach to low velocity impact analysis is to develop a 
specialty finite element program [4-14] that can account for the nonlinear elastic 
force— displacement behavior in the contact region. Often a triangular of half— sine 
waveform is assumed [1, 4] to represent the transient impact force history. In the work 
presented here, a more straightforward method of impact modeling is presented, using the 
general purpose structural analysis program NASTRAN. A nonlinear elastic contact law is 
used to model the local contact behavior, and the resulting impact force is calculated based 
on the relative displacement between impactor and target. 


BACKGROUND 

Structural damage due to impact invariably initiates in the immediate vicinity 
of the impact. Therefore, it is important that the local stress field in the region of contact 
be calculated accurately. The Hertzian contact law [15] is as an elasticity-based 
force— displacement relationship that describes contact between two elastic bodies. The 
Hertzian contact law is given by: 


F = Ka" 


( 1 ) 


where 

F = elastic contact force 
K = contact stiffness 
n = exponent 

and 

a = relative displacement (indentation) between impactor and target 
= Ui — u t (i = impactor, t = target) 

The exponent n was shown in reference [15] to have the value of 3 / 2 - In this application, 
force, displacement and indentation change with time. 
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During low velocity impact, structural damage to the target is negligible, and 
areas of the structure remote from the impact deform in a linear elastic manner. An 
efficient finite element model, therefore, would combine a linear elastic model of the global 
structure with the prescribed non-linear elastic behavior at the point of contact with the 
projectile. The nonlinear force— displacement relationship in equation (1) is incorporated 
into a linear elastic finite element model (MSC/NASTRAN transient solution 27, 
COSMIC/NASTRAN transient solution 9) by using a NASTRAN transfer function 
definition and nonlinear analysis capability. In the following section, Hertz’ contact law is 
discussed. The next section describes a straightforward method of incorporating this 
contact law into NASTRAN. Impact loading of two different structures is then analyzed. 
The first problem considered is a one— dimensional rod of uniform cross section impacted 
longitudinally. The second is an orthotropic plate under transverse impact. Calculated 
results are compared with experimental test data for both cases. 


CONTACT LAW 

In reference [16], the derivation of the Hertzian force— displacement relationship 
is given for two spherical isotropic elastic bodies of radius r! and r 2 in contact. The 
resulting equation is given by: 

F = K a 3 / 2 (2) 

where 


is the contact stiffness and 




ri r 2 

CN 

ri + r 2 

< 
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( 3 ) 


( 4 ) 


where E^. and v- are the elastic modulus and Poisson ratio, respectively, and the subscripts 

1 and 2 refer to each of the spheres. When a spherical impactor contacts a flat target, 
equation (3) simplifies to 


K 




hi 

hi + 



( 5 ) 


where i and t represent the impactor and target, respectively, and hj and h t are given by: 
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( 6 ) 


*i = r 

1 - Vi 

kt (7) 

1 -I/ t 

In equation (2), a is the local indentation at the contact point, shown 
schematically in figure 1. We have: 


a = uj-u t (8) 

where a is the relative local displacement between impactor and target at the point of 
contact. 


NASTRAN Implementation 

NASTRAN allows for easy implementation of the contact law and is readily 
available to many users. This makes NASTRAN a excellent finite element program for 
incorporating this method of impact modeling. A unique feature of the this impact 
analysis is that the impact force is not a pre— defined input parameter, but rather is 
calculated during the analysis, based on the relative displacement between impactor and 
target, and using the prescribed force— displacement relationship. This methodology allows 
for the unloading of the structure when the target and the impactor are no longer in 
contact and for secondary impacts in the event that they would again come into contact. 
The contact law was incorporated into the NASTRAN finite element model as follows: 

The impactor is modeled as a lumped mass just touching the target at time t=0 
and with an initial velocity toward the target. The difference between the displacement of 
this lumped mass and the displacement of the target is the indentation, a. The contact 
law is prescribed with the transfer function (TF) card, and the nonlinear force (NOLIN3) 
card. The TF card acts as a dynamic multipoint constraint, relating the displacement, 
velocity and acceleration of several independent degrees of freedom to a single dependent 
degree of freedom. In this case, only displacement relationships were needed. The 
dependant degree of freedom is calculated on the TF card as follows: [17]. 


ii 

(B 0 + B t p + Bjp^p + £(A J 0+ Aj P + A J 2 p2)u J ind = 0 


j 


(9) 


where 


B 0 , Bj, B 2 = the coefficients for the dependent degree of freedom 
Aq, Aj, A J a = the coefficients for the independent degrees of freedom 
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u dep = displacement of the dependent degree of freedom 

u ind = displacements of the independent degrees of freedom 
n = the number of independent degrees of freedom 
p = the differential operator, ^ ; and p 2 = — ^ 


For this analysis, the equation would appear: 

(l.O)Uextra point "h ^( _ 1-O)uimpactor + (l.O)utargetj = 0 


( 10 ) 


that is 


n = 2 


®i> Bj, Aj, Aj 

B ft 


0.0 (j = l,n) 
1.0 


Aj = -1.0 

A l = 1.0 


The resulting equation defines the dependent degree of freedom as the indentation at each 
time step of the transient analysis. The value of the indentation is assigned to an 
EPOINT. The EPOINT, or extra point, is a nonstructural variable that is used to store 
the value of the indentation. The EPOINT is provided as input to the NOLIN3 card. 

The NOLIN3 card is the means of applying the time— dependent nonlinear load 
based on the indentation. The NOLIN3 card has the form: 


where 


P(t) 


S[x(t)] A , x(t) > 0 
0 , x(t) < 0 


( 11 ) 


P(t) = the resulting nonlinear force 
S = a scale factor 

x(t) = the displacement or velocity of a degree of freedom 
A = an amplification factor 
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In modeling of the impact, we define x(t) to be the displacement of the EPOINT, S to be 
the Hertzian stiffness, and A to be 3 / 2 , as given in equation (2). Recall that the 
displacement of the EPOINT is really the indentation as obtained from the TF card. The 
resulting function then has the form: 


F(t) 


K[a(t)] 3 /’, o(t) > 0 
0 , o(t) < 0 


( 12 ) 


Note that when a is less than or equal to zero (target and the impactor are no longer in 
contact), the contact force is also zero. Two NOLIN3 cards are used, one to apply the 
impact force to the target and the other to apply the same force to the impactor in the 
opposite direction of its initial velocity. This methodology allows the impactor to slow 
with increasing impact force and eventually to unload the target as the impactor begins to 
travel in the opposite direction, away from the target. Since the forcing function is not an 
input parameter,, changing the properties of the target or the impactor does not require 
assuming a new force-time relationship but only the recalculation of the contact stiffness, 
K, and changing the appropriate material property cards. 


RESULTS 

One Dimensional Rod 

The first problem analyzed is the longitudinal impact of a steel ball on a long 
aluminum rod of constant cross section. Geometry and material properties of the impactor 
and target are given in figure 1. The problem was modeled using 144 l-<limensional rod 
elements with each grid point having a single longitudinal degree of freedom. Two 
additional degrees of freedom were used to model impactor and the extra point, resulting in 
a total of 147 degrees of freedom. A single lumped mass with an initial velocity was used 
to represent the impactor. The Hertzian force— displacement relationship in equation (1) 
was prescribed using the NASTRAN NOLIN3 card, as shown in the example input deck in 
the appendix. 

The impact force history obtained from the analysis compares well with 
experimentally determined values [18], as shown in figure 2. The calculated strain response 
at the midpoint of the target bar is compared with measured values in figure 3. The sign 
reversal of the second pulse is caused by the reflected tensile stress wave generated by the 
incident compressive wave reaching the free end of the bar [19]. 

Some insight into the timing and the location of the impact-induced structural 
failure can be gained by tracking the distribution of energy in the impactor and the target, 
as shown in figure 4. The energy balance can be expressed as: 

Utot = KEi + SEj + KE t + SE t (13) 


where 


6 



Utot = total energy in system 
KEi = impactor kinetic energy = -y- mvj 

SEi = impactor strain energy = J F(a) da = — g— Ka^ 2 
KEt = target kinetic energy 

n 

■s+ 

SEt = strain energy of target 
1 2 

— y kj ( — £j +1 ) (n = number of elements) 



j — (n = number of elements) 


(14) 

(15) 

(16) 


(17) 


The total energy in the system, Utot, is divided between the kinetic energy and 
the strain energy of the target and the impactor in a time— varying manner. Because 
damping effects are not considered, the total system energy is constant and equal to the 
initial kinetic energy of the impactor. The strain energy of the impactor is non-zero only 
during the contact interval (0 < tc/L < 0.4, where t is time, L is the length of the bar, and 
c is the wave speed in the bar) and peaks when the contact force is greatest, approximately 
halfway through the contact interval. The kinetic energy of the impactor decreases rapidly 
as the impactor slows during contact with the target. Eventually, at tc/L £ 0.25, the 
impactor velocity (and therefore its kinetic energy) decreases to zero and "the elastic 
rebound begins. The kinetic energy of the impactor never returns to its initial level, 
however, because approximately 80% of the energy has been transferred to the target in the 
form of strain energy and kinetic energy. The strain and kinetic energies in the target both 
increase rapidly during the contact with the impactor and remain constant after contact 
has ended (tc/L > 0.4). Both strain and kinetic energies maintain equal and constant 
values until the compressive stress wave generated by the impact reaches the far end of the 
free— free bar (tc/L = 1.0). A tensile stress wave is generated when the compressive pulse 
reflects from the stress free boundary [19]. The superposition of the incident and reflected 
pulses momentarily leaves the bar stress— free, and the strain energy decreases to zero. The 
kinetic energy simultaneously increases, maintaining a conservation of total energy. The 
reflection process is repeated at tc/L = 2.0, when the reflected pulse returns to the 
impacted end of the bar. Similar energy dissipation diagrams may prove useful in 
analyzing dynamic failure of more complex structures. By tracking the distribution of 
energy and assuming a strain energy based failure criterion, it maybe possible to locate 
areas of failure in the target and predict the time of failure of these areas. 


Composite Plate 

The low velocity transverse impact of a composite plate made from Scotchply 
1003 prepreg [20] is now analyzed. The problem is depicted schematically in figure 5, and 
is described in detail in references [9, 10j. A modified Hertzian contact stiffness has been 
proposed [5] for application to composite materials. Specifically, equation (7) is replaced 
by 


kt = E33 


(18) 
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where E 33 is the transverse (out-of-plane) modulus of the composite plate. Plate 
membrane and bending stiffness material properties were calculated using the COBSTRAN 
(Composite Blade Structural Analyzer) computer code [21] which calculates elastic moduli 
of composite materials from known constituent properties and laminate ply orientations. 

A uniform square mesh of QUAD4 elements was used to model the 15.24 cm * 
15.24 cm (6 in * 6 in) target plate. A mesh convergence study was performed to establish 
the degree of mesh refinement necessary to arrive at a numerically converged solution. 
Three different meshes were considered, 25 * 25, 49 * 49, and 61 * 61 elements. Of these, 
the latter two produced essentially the same strain response for a given impact velocity and 
were therefore considered to be converged solutions. The results presented here were 
therefore calculated using the 49 * 49 element model. Five degrees of freedom (u x , u y , u z , 
0* and 0 y ) were used at each nodal point, giving the model a total of 11510 degrees of 
freedom. The problem was solved on a Cray XMP in 52 CPU minutes. 

The impactor used in the tests [9, 10] was a uniform 2.54 cm (1 in) long, 
blunt-ended steel rod of radius 0.048 cm ( 3 / 16 in). In the analysis, a contact radius of 
0.047625 cm ( 3 /u in) was assumed in the Hertzian contact stiffness calculations. The 
calculated impact force history is shown in figure 6. Although no direct measurement of 
the impact force was obtained experimentally, the contact time was measured [8] and fo un d 
to be 204 microseconds. This is in good agreement with the calculated result. Figure 6 
also shows that a secondary impact occurs during the latter half of the contact interval (t 
= 175 pse c), probably due to the vibration of the target plate during contact with the 
impactor. 

The resulting displacement response of the plate is shown in figure 7, where it 
has been assumed that no damage occurs in the target during contact with the impactor. 
This assumption is valid based on the available test data. Ultrasonic C— scans of the 
specimens after impact indicate that this level of impactor kinetic energy (10 Joules) is 
very near the threshold energy level required to cause damage [10] in specimens of this 
layup. As a result, very little damage occurs at this impact velocity. The anisotropic 
bending stiffness of the target is evident from the elliptical displacement contours, as the 
flexural disturbance travels faster in the stiffer direction, as shown in figure 7. 

The strain response at gage A is compared to the calculated response in figure 8. 
The two curves are similar in amplitude and duration but the calculated strain appears to 
lag the measured values by approximately 25 microseconds. This may be due to the 
difficulty in establishing experimentally the precise time at which contact occurs based on 
strain gage readings taken at some distance from the point of contact. The comparison 
shown in figure 9 for gage B likewise shows a time shift of approximately 25 microseconds 
between the measured and the calculated response. The amplitude and duration of the 
calculated strain response correlate quite well with the measured signal. 
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SUMMARY 


A simple means of modeling low velocity, non— damaging impact using 
NASTRAN was demonstrated. A nonlinear elastic contact model was included in the finite 
element analysis using NASTRAN transfer function definitions and nonlinear analysis 
capabilities. The same contact law was used to define the force-indentation relationship 
for two different impactor/target combinations. Results in both cases showed that the 
impact force and resulting transient structural response of the target compared well with 
experimentally measured values. Future work will include the effects of damage and the 
resulting progressive loss of structural stiffness that occurs during higher velocity impact 
events.- 
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DAN TROWBRIDGE 


ANALEX - STRUCTRAl MECHANICS BRANCH 


ID TRANS. LOAD 
APP DISP 
TIME 60 
SOL 9 
CEND 

TITLE - COSMIC: TRANSIENT RESPONSE ANALYSIS: HERTZIAN IMPACT FF 
SUBTITLE - 36 H AL. ROD 5/8 STEEL BALL V0=62.1 IN/S 

LABEL * ROD - * < IMPACT 

$ NONLINEAR LOAD 

NONLINEAR - 5 

$ * INITIAL CONDITIONS SET 

IC * 1 

TFL*111 

SPC - 4 

TSTEP - 7 

$ OUTPUT STUFF 

SET 30 - 1.72.73.999.1001 
NLLOAD * 30 
STRESS(PRINT) = 30 
DISP (PR I NT) * 30 
BEGIN BULK 

$ 

$ EXTRA POINT - INDENTATION 

EPOINT. 1001 

GRID. 999. ,-0.3125,0.0,0.0 
GRID, 1 , .0.0. 0.0, 0.0 
*(144), *(1),*. *(0.25),= 

CROC , 1 , 1 , 1 , 2 
«<143)..(1).«..(1),*(1) 

$ LUMP MASS OF IMP AC TOR 

CONM2 . 200 .999,0,9. 587-5 ,0.0,0. 0,0.0,, +CON2-2 
+CON2-2.3.745-6, ,3.745-6. , ,3.745-6 
$ MATERIAL PROPERTIES 

PROD , 1 ,11,0.196,6. 14-3,0.25 

MAT 1,11,10. 0-f 6 , ,0.33,2. 5-4 +MAT1-1 

+MAT1-1 , 35 . 0E6 . 36 . 0E6 , 27 . 0E6 
$ BOUNDRY CONDITIONS 

SPC1. 4. 23456,1, THRU, 145 

$ REMOVE DEGREES OF FREEDOM FROM IMPACTOR 

SPC1 , 4, 23456, 999 

$ TRANSFER FUNCTION TO DEFINE INDENTATION 

TF. Ill ,1001 ,0.41 .0,0. 0,0.0,, ,+TF-l 

+TF— 1 ,999,1 ,-1 .0,0. 0,0.0, ,, .+TF-2 

+TF-2. 1 .1,1. 0,0. 0,0.0 

$ TIMING 

TSTEP. 7, 2500, 2. 0-7, 25 

$ LOAD DEPENDENT ON DISPLACEMENT OF IMPACTOR 

NOLIN3 , 5, 1. 1, 6.2446. 1001, 1. 1.5 
$ SLOW DOWN IMPACTOR 

NOLIN3, 5, 999, 1, -6.24+6. 1001. 1, 1.5 

$ INITIAL CONDITIONS: IMPACTOR VELOCITY - 62.1 IN/SEC 

TIC, 1 .999.1 ,0.0,62.1 

ENDDATA 
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DAN TROWBRIDGE 


ANALEX - STRUCTRAl MECHANICS BRANCH 


ID IMPACT. PLATE 
APP OISP 
TIME 120 
SOL 27 
CEND 

TITLE - IMPACT OF PLATE 49X49 : CENTERED ELEMENT 
SUBTITLE - TRANSIENT ANALYSIS: FIXED-FIXED: NO SYMMETRY 
SPC - 1 
IC - 3 

NONLINEAR - S 
TSTEP - 1 
TFL - 111 

SET 15 - 999.2525.2526.2625.2626 
NLLOAD - 15 

SET 20 - 2525.3525.4125 
STRESS - 20 
BEGIN BULK 

$ • ••• EXTRA POINT TO HOLD INDENTATION ................. 

EPOINT. 10001 

$ .... IMPACTOR 3/8 IN DIAMETER 

GRID. 999. .0.0.0.0,-0.1875 

C0NM2. 200. 999. 0.8. 096-5. 0.0. 0.0. 0.0, . +CON2-2 
+C0N2— 2 . 7 . 459—6 . .7.459-6, , . 1 .423-6 


$ •••••. GRIDS AND CQUAD4 ELEMENTS DEFINING THE PLATE GO HERE ... 

$ MATERIAL PROPERTIES. . . MAT2 CARDS GENERATED BY COBSTRAN 


PSHELL.1 .101.0.15.201 .1 .0 

MAT2.101 , 4.3E+06 , 2 . 9E+05 .— 1 . 7E— 03 , 2 . 8E+06 .-3 . 4E-02 . 5 . 7E+05 . 1 .8E— 04.+A1 01 
+A101 .5.8E— 06.8.9E-06.5.0E-13 

MAT2.201 . 5 . 7E+06 , 2 . 9E+05 , -1 .9E-04, 1 . 4E+06 ,-3 . 8E-03 . 5 . 7E+05 
S BOUNDRY CONDITIONS 

SPC1 , 1. 123456. 101, THRU. 150 
SPC1 . 1. 123456. 5001. THRU, 505e 
SPC1 . 1. 123456, 101 
-.-.-.•100 
-48 

SPC1. 1. 123456, 150 
-.-.-.•100 
=48 

SPC1 . 1. 12456. 999 
GRDSET 6 

$ TIME STEP INFO 

TSTEP. 1.2000. 1.0-7. 10 

S LOAD DEPENDENT ON RELATIVE DISPLACEMENT OF IMPACTOR 


N0LIN3, 5.2525. 3. +1.945+5 

. 10001, 

0, 

1 *5 

N0LIN3. 5.2526, 3. +1.945+5 

, 10001. 

0, 

1 .5 

N0LIN3. 5.2625. 3. +1.945+5 

. 10001. 

0. 

1.5 

N0LIN3. 5.2626. 3. +1.945+5 

. 10001, 

0. 

1.5 

% SLOW DOWN 

IMPACTOR 



N0LIN3. 5. 999. 3.-7.779+5 

, 10001, 

0. 

1 .5 


S TRANSFER FUNCTION TO CALCULATE INDENTATION 

TF.111 .10001 .0.+1 .0,00.0,00.0, . ,+TF-l 

+TF-1 .999.3.-1 .0.00.0,00.0, . . , +TF— 2 

+TF-2. 2525. 3. +0.25. 00. 0,00.0. . . .+TF-3 

+TF-3. 2526, 3. +0.25, 00. 0,00.0, , . .+TF-4 

+TF-4. 2625. 3. +0.25, 00. 0,00.0, . . .+TF-5 

+TF-5. 2626. 3. +0.25. 00. 0,00.0 

* INITIAL CONDITIONS: IMPACTOR VELOCITY - 1470 IN/SEC (122.5 FT/SEC) 

TIC. 3. 999. 3. 0.0, 1470.0 

ENOOATA 
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STRAIN, p h/h 


1.59 CM (0.S25 in.) 
DIAMETER V Q 
STEEL BALL O 


91.44 cm (36 IN.) 



' V. '-53 

1.27 CM (0.5 in.) DIAf€TER ALUMINUM ROD 


V 0 = 158 CM/S (82.1 IN./S) 





BAR. 


STRAIN ENERGY (BAR) 

— — — KINETIC ENERGY (BAR) 

STRAIN ENERGY (BALL) 

1.2 — KINETIC ENERGY (BALL) 

TOTAL ENERGY 



0 .5 1.0 1.5 2.0 2.5 

tc 


L 

FIGURE 4. - ENERGY DISTRIBUTION FOR LONGITUDINAL BAR 
IMPACT PROBLEM. 
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STRAIN. V*/H 



FIGURE 5. - COMPOSITE PLATE IMPACT SPECIMEN CONFIGURATION 
(REF. 9. 10). 



FIGURE 6. - CALCULATED FORCE HISTORY FOR TRANSVERSE 
IMPACT OF COMPOSITE PLATE. 




60 msec 



200 msec 


FIGURE 7. - CALCULATED DISPLACEMENT RESPONSE FOR COMPOSITE 
PLATE. 
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